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We derive analytical solutions for the system of two ultracold spin-polarized fermions interacting 
in p wave and confined in an axially symmetric harmonic trap. To this end we utilize p-wave pseu- 
dopotential with an energy-dependent scattering volume. This allows to describe the scattering in 
tight trapping potentials in the presence of scattering resonances. We verify predictions of the pseu- 
dopotential treatment for some model interaction potential, obtaining an excellent agreement with 
exact energy levels. Then we turn to the experimentally relevant case of neutral atom interactions 
in the vicinity of a p-wave Feshbach resonance. In the framework of the multichannel quantum- 
defect theory we derive relatively simple formula for an energy-dependent scattering volume, and 
later we apply it to investigate the energy spectrum of trapped atoms close to the p-wave Feshbach 
resonance. 
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I. INTRODUCTION 

In the regime of ultracold temperatures atomic col- 
lisions are typically dominated by s-wave scattering. 
The situation changes in the presence of scattering reso- 
nances, which can strongly enhance the contribution from 
higher partial waves. This has been demonstrated in ex- 
periments by employing Feshbach resonances to enhance 
p-wave interactions in Fermi gases [H, H, H Q and Fermi- 
Bose mixtures Q ■ So far realization of p-wave superfiuids 
is hindered by the presence of strong collisional losses in 
the vicinity of resonances, and, as a consequence, short 
lifetime of the ultracold gas. Similar losses affect the 
lifetime of p-wave Feshbach molecules, that has been re- 
cently created [f| 0, HI and theoretically investigated @ . 
Realization of p-wave superfiuids is of particular concern 
because of their very interesting properties, e.g. rich 
phase diagram with a variet y o f classical and quantum 
phase transitions 0, EH, E 111, • 

Theoretical description of ultracold gases typically as- 
sumes modeling of an atom-atom interaction in terms 
of a contact potential. In the case of s-wave interac- 
tions this is typically done by means of Fermi pseudopo- 
tential [HI, Analogous pseudopotential for higher 

partial waves has been developed several years ago by 
Huang and Yang [l(| [TtJ , however their original deriva- 
tion turned out to be incorrect [H, and resulting 
pseudopotential valid only for s-wave scattering. Dis- 
covery of an error in the derivation of Huang and Yang 
has been a motivation for a series of papers propos- 
ing several forms of p seudopotentials for higher partial 
waves [H [H [H M, [H, El Si, HI. Such pseudopo- 
tentials have been applied, for instance, to confirm the 
existence of geometric resonances [2(| H3] in the scatter- 
ing of fermions in quasi-two-dimensional systems [23| or 
to investigate trap-induced shape resonances of p-wave 
interacting atoms [24j . 

In this paper we apply the pseudopotential method 
to develop analytical solutions for two spin-polarized 
fermions interacting in p-wave and confined in the axi- 
ally symmetric harmonic trap. Some of the results for 



pancake-shape traps have been already presented in our 
earlier work [23j |. Here we present much more detailed 
analysis, including the analytical formulas for traps of 
arbitrary aspect ratio. We verify the accuracy of the 
pseudopotential treatment by comparison with exact en- 
ergy levels for a model square-well potential. Finally we 
consider a relevant case of atom-atom interactions tuned 
by the p-wave Feshbach resonances. Our results are of di- 
rect interest for studies of spin-polarized Fermi gases in 
optical lattices Q, where tightly confined pairs of atoms 
can be obtained and p-wave Feshbach molecules can be 
created by applying time-dependent magnetic fields, sim- 
ilarly to s-wave Feshbach resonances [H [2lJ |3(| ■ 

An important part of our analysis constitute rigorous 
derivation of pseudopotentials for higher partial waves. 
Here, for p-wave interactions we utilize slightly modified 
version of zero-range potential than introduced in [23| . 
which is more convenient for our analytical calculations 
in the trap, and allows for simple generalization for par- 
tial waves I > 1. Because of the variety of pseudopo- 
tentials in the literature, we devote a small section on 
comparison of our pseudopotential to other formulations, 
discussing the conditions when they become equivalent. 

In contrast to s-wave interaction where the scattering 
amplitude acquires a constant value at low energy, and 
can be characterized with a single parameter (e.g. s- 
wave scattering length), pseudopotentials in higher par- 
tial requires inclusion of an energy-dependent amplitude 
for modeling of resonant interactions. Existence of the 
centrifugal barrier makes the higher partial wave reso- 
nances intrinsically narrow and sensitive to the energy 
of the colliding particles. In our approach we employ 
multichannel- quan tum defect theory (MQDT) for ultra- 
cold collisions [3lL I32I [33[ , that allows for an exact treat- 
ment of two-body scattering in the presence resonances 
and yields relatively simple expression for the energy- 
dependent amplitude of the pseudopotential due to the 
p-wave Feshbach resonance. 

The structure of the paper is as follows. In section 
III Al we present detailed derivation of the pseudopoten- 
tial for all partial wave interactions. Section IIIBI con- 
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tains a comparison of different pseudopotentials proposed 
in the literature. Section III CI presents simpler form of 
pseudopotential, that can be used when interaction takes 
place only in a single partial wave. Reader not interested 
in details of derivation, may skip sections III AIIII Cl and 
jump directly to section Til Dl where we introduce p-wave 
pseudopotential used throughout this paper, and also, 
for the sake of illustration, a pseudopotential for d-wave 
interactions. Section HTT1 is devoted to analysis of two in- 
teracting spin-polarized fermions in a harmonic trap. In 
Section HITS we derive some general analytical result for 
the energy spectrum. Specifying the form of the energy- 
dependent scattering volume, we apply these results to 
investigate the energy spectrum for a model square-well 
interaction (Section lHI Bp and for an atom-atom interac- 
tion tuned with a magnetic Feshbach resonance (Section 
MI C|) . We summarize in Section lTVl presenting some con- 
clusions. Finally Appendices [AllCl present some technical 
details related to derivation of pseudopotentials. 



II. PSEUDOPOTENTIALS FOR HIGHER 
PARTIAL WAVE INTERACTIONS 



pseudopotential comes from ni(kr) that is singular at 
small r: m(x) = -(21 - l)\\/x l+1 + ©(l/a;^ 1 ). We have 
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For terms of the order of and higher in the ex- 

pansion of ni(kr), we observe that the radial part of the 
Laplace operator produces terms that are exactly can- 
celed by terms obtained from the multiplication by k 2 . 
On the other hand, the most singular term l/r l+1 gives 
no contribution at r ^ 0, and the only nontrivial part 
comes from the radial part of the Laplace operator ap- 
plied to l/r l+1 at r = 0. The latter can be calculated on 
the ground of the theory of generalized functions. 

In our derivation we adopt the Hadamard finite part 
regularization of the singular function l/r k , which can 
be defined as M, M 



A. Corrected pseudopotential of Huang and Yang 

In this section we present detailed derivation of the 
generalized pseudo pot ential for all partial waves, intro- 
duced in Refs. (l9l |23||. It can be obtained in a manner 
similar to the original derivation of Huang and Yang [l6| , 
where the only difference is the rigorous treatment of the 
delta-function derivatives within the distribution theory. 
We start from the Schrodinger equation in the relative 
coordinate 
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(A + /j 2 )*(r) = V(r)#(r), 



(1) 



where k 2 — 2/j.E/h 2 , and \i denotes the reduced mass. 
We assume that potential V(r) is central and has a fi- 
nite range [56[ . Outside the range of potential, the wave 
function exhibits the following asymptotic behavior 



^( r )=E E R i(r)Yim(i), 



(2) 
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Here, Yi m (r) are spherical harmonics, the radial wave 
functions Ri(r) = Ci m (ji(kr) — tan£;n/(fcr)) are linear 
combinations of spherical Bessel and von Neumann func- 
tions ji(r) and ni(r), respectively, Ci m are coefficients 
that depend on the boundary conditions at r — > oo, and 
the phase shifts Si are determined by the potential V(r). 
Our goal is to replace V(r) by zero-range potential V pe (r), 
which acts only at r = and gives the same asymptotic 
function ^ a (r), as the real potential V(r). Following 
Huang and Yang [l6[ , we investigate the behavior of the 
l.h.s. of Eq. flU at r — > 0, substituting * = * Q . Since 
ji(kr) is regular at r = 0, the only contribution to the 
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for k integer and positive. This definition follows simply 
from substraction of diverging terms in the integral of the 



type J dxcj)(x)/x k 



(see Appendix A for details) . 



The presence of the Heaviside function H(r) is related 
to the fact that the radial wave function is defined for 
r > 0. We stress that regularization procedure is not 
unique (see for instance [3J]), however all regularizations 
must give the same results when integrated with a test 
function behaving like r k for small r. Regularization Q 
has the following properties [H| 
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where d/(dr) represents generalized, distributional deri- 
vative, and S^(r) denotes the n-th derivative od the 
delta function. When the distribution 6^ is applied to 
a test function <j)(r) that vanishes at r — > as r a with 
a > 1- it yields 

(<p,6W(r)) = -kU^^-^M/r), lim ^ < oo, (7) 

r— >o r 

which follows easily from integration by parts. Combin- 
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ing ([6]) and $7} we obtain 
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(8) 

By applying ([7]) twice we have extracted 1/r 2 term, 
which is convenient when calculating matrix elements in 
spherical coordinates, because of the additional r 2 fac- 
tor arising from the Jacobian. Finally, we express the 
coefficients Ci m in terms of the radial wave function: 
C lm = ((2l + l)\\/[k\2l + iy]{d 2l + 1 r l + 1 R l (r)) r=Q . The 
regularization operator d 2l+1 r l+1 removes the singular 
contribution from n;(Ar), whereas the prefactor is de- 
termined by the small argument behavior of the ji(kr): 
ji(x) = x l /(2l + 1)!! + 0(x l+2 ). This leads to the pseu- 
dopotential [23j 
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where the angular integral over f2' acts a projection op- 
erator on a state with a given quantum number I, and we 
have performed the summation over m resulting in the 
Legendre polynomial p. Here, gi = h 2 ai(k) 2l+1 /(2/x), 
and ai(k) denotes the energy-dependent scattering length 
for partial wave I: ai{k) 2l+1 = - tan 5i(k)/k 2l+1 jHEEzI- 
We observe that I = component yields the familiar 
s-wave contact potential 2irh 2 ao(k) / [iS(r)d r r, whereas 
components with I > differ from the result of Huang 
and Yang by a prefactor (21 + I) /(I + 1) [l7|. 

Finally we note that for I > 1 the Bessel function ni(x), 
apart from the leading order singularity l/x l+1 , contains 
also lower order singularities, which assuming our regu- 
larization, in principle give rise to the lower-order deriva- 
tives of the delta distribution for Z-th component of the 
pseudopotential. Such terms, however, do not give any 
contribution for functions behaving like r l for small r, 
and they do not appear when the other regularizations 
[25[ are applied. 



B. Other representations of pseudopotential 

Now we show, the equivalence of ([9]) to other zero- 
range potentials derived recently in the literature (nil . l2(il . 
[25j . To this end we evaluate matrix elements of V ps , mak- 
ing some assumption about the form of the wave function 
$(r) acting to the l.h.s. of Q. First, assuming that the 
l-w&ve components 4>i(r) of <&(r) behave as 4>i(r) ~ r a 
with a > I at small r, we can apply ||7J) recursively, to 
obtain 



((MW(r)) = (-l) l ll(^,6(r)/r l ), lim^(r) 

r— >0 
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This combined with the identity 6(r) = 8(r) / (Airr 2 ) re- 
sults in the pseudopotential 1!) 
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The delta-shell pseudopotential of Stock et al. |2(j, dif- 
fers from the latter representation by use of the delta- 
shell distribution lim s _>o 8(r — s)/s 2 , instead of S(r)/r 2 . 

Another useful representation of the pseudopotential 
was introduced by Stampfer and Wagner [25| 
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where the partial differential operator Y\ m (d) is obtained 
from the harmonic polynomial r Yj m (r), by replacing the 
Cartesian coordinates Xk with the partial derivatives d Xk . 
The operator Yi m (d), introduced by Maxwell [38|, has the 
following properties 

Y lm (d)r x = A(A — 2) ... (A — 2/ + 2)r A -%„(f ), (13) 

= (21+1 )11, 
r=o Air 
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The equivalence of pseudopotentials 
from the following identity 



and (Hill follows 
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which we prove in the Appendix B. 



C. Pseudopotentials for interactions in a single 
partial wave 

Quite often the atom-atom interaction occurs domi- 
nantly in a single partial wave, e.g. in the vicinity of 
resonances. In such cases we can consider only a single 
partial wave, and derive simpler form of the pseudopo- 
tential. To this end, we express the projection operator 
on a partial wave I, that appears int the square brackets 
in Eq. (TT2")) , in terms of the differential operator Yi m (d). 
We observe that the only non vanishing contribution to 
a2'+V' +1 /dfi'Y4(r / )#(r / )l , comes from the terms 

J r'=0 

of the symmetry r Y/ m (f'). Therefore we can rewrite the 
regularization part of the pseudopotential, introducing 
an additional differentiation and multiplication by r', 
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In a similar fashion we can transform the projection op- 
erator 



1 d l 



which does not affect the part of the wave function pro- 
portional to (r') Yi m (r'), as can be easily verified apply- 
ing the property (|15[) . This yields the following represen- 
tation of the zero-range potential for partial wave I 

VKr)tt(r) = Z ■ J 9lYi m (d)6(r) 
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(21 + 1)1 
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r' 2 ' + V^(a')*(r') 
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Pseudopotential (|18p involves only differential operators, 
that are typically much easier to handle than the angular 
integration projecting on a given partial wave I in (|12[) . 
We note that representation (|18[) is equivalent to (TT2"|) 
only if the interactions occur in a single partial wave I, 
and the resulting wave function exhibits 1/r or r' be- 
havior at small distances. 

We observe that the summation over quantum number 
to involves only the differential operators Y/ m (<9). The 
sum over to can be performed using an identity that di- 
rectly follows from the summation formula for the spher- 
ical harmonics Yi m (r) (see Appendix C for derivation) 

£ Y lrn {d)Y? m {d') = 

m=—l 

= £c fc (V- V')^ 2fc V 2fe (V') 2fe . (19) 

k=0 

Here, are coefficients of the Legendre polynomial 

m 

Pi(x) — J2k~o c kX l ~ 2k , that are given by 
(-l) k (2l - 2k)\ 



Ck = 



2 l k\{l - k)\{l - 2k)V 
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and [a;] denotes the largest integer, smaller or equal to x. 
In the following section we will apply representation (Tig]) 
and summation formula (|19| to derive relatively simple 
expression for the p- and d-wave zero-range potentials. 



D. Pseudopotentials for p and d wave in Cartesian 
coordinates 



In the problems involving interacting atoms, very often 
an external potential does not posses a spherical symme- 
try, and the total angular momentum is not conserved. In 
such systems modeling of interactions in terms of pseu- 
dopotentials ©, (fTTj) . or (fT2"|) , containing explicit inte- 
gration with spherical harmonics, is not convenient. One 



of examples is the exact treatment of few-body systems 
confined in tight anisotropic harmonic trap, that is rel- 
evant in context of experiments on ultracold atoms in 
optical lattices. In such cases, the pseudopotential l(T5j) 
is far more easier to apply. 

In this section we discuss pseudopotentials for p- and d- 
wave interactions, writing down explicitly the differential 
operators in the Cartesian coordinate system. Applying 
the summation formula (fT9]) to the p-wave part of fTS"]) 
we arrive at 



V p (r) = p - V <Kr)^r 3 V , 

fx or 6 



(21) 



where the symbol V (V) denotes the gradient opera- 
tor that acts to the left (right) of the pseudopotential, 
and a p is the p-wave energy-dependent scattering length: 
a p (k) 3 — — tan(5i(fc)/fc 3 . This form of the pseudopoten- 
tial with a scalar product of the two gradient operators 
have been first introduced by Omont [3^|, and later by 
Kanjilal and Blume [2l|, however, in the former the reg- 
ularization operator was missing, whereas in the latter 
the regularization was different than ours. 

Similar procedure applied to c?-wave part leads to 



V d (r) = 
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where D 
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^SijSki and ad is the d-wave scat- 
tering length: a^ik) 5 = — tan(S 2 (fc)/fc 5 . The zero-range 
potentials (|2"Tj) and (f2"2"| presented here are slightly dif- 
ferent from those introduced in Ref . [23| , differing in the 
form of the regularization potential and the position of 
the differential operator acting to the right. One can 
easily show that they are equivalent, whereas the present 
derivation is much easier to generalize to arbitrary par- 
tial wave, resulting in a formula (|19| that can be used to 
generate pseudopotentials for I > 2. 

Finally, we note that s-wave component of the pseu- 
dopotential in the representation (|18[) contains only the 
regularization operator d r r, hence it does not require fur- 
ther transformations, and it can be directly applied in the 
problems with anisotropic external potentials We 
stress that the use of s-wave zero-range potential with- 
out the projection operator is valid as long as interactions 
takes place only in s wave, and the radial wave function 
can exhibit only 1/r singular behavior at small r. 



III. TWO ATOMS INTERACTING IN p-WAVE 
AND CONFINED IN A HARMONIC TRAP 
POTENTIAL 

A. Analytical results for the zero-range potential 

In the following we consider two spin-polarized fer- 
mions confined in an anisotropic harmonic trap. At low 
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energies the interaction is dominated by the p-wave, and 
we model it using the zero-range potential (j2"Tj) . We cal- 
culate the energy spectrum of the system, that can be 
obtained solving the Lippmann-Schwinger equation for 
the wave function \&(r) of the relative motion 



tf(r)= /dV 'G(r,r')V p (r')#(r'), 



(23) 



that for the discrete spectrum does not contain the source 
term. Here, Ge(t,t') = (r\(E — i?) _1 |r') is the single- 
particle Green function, and H is the Hamiltonian in- 
cluding the external potential. For anisotropic harmonic 
potential Ge(r, r') can be represented by the following 
integral [57j ] 
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which is convergent for energies smaller than the energy 
of zero-point oscillations Eq = h(ui x +u)y+ui z )/2. For E > 
Eq the Green function can be determined by analytic 
continuation of ([24l) . Here, d = y/h/Jpuj), % = oJk/uj, 
u> denotes some reference frequency, and product in (|2^|) 
runs over k = x,y, z. Inserting the pscudopotential (|2ip 
into (|2"3")) . we obtain the following set of equations 



C k = p -Y,Au{E)C u 



where matrix elements Am(E) are given by 



d 3 r 3 d 2 Ce{v r / } 

dr 3 dxkdx't 
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and coefficients Ci = [d 3 r 2 d Xl ^(r)^ r Q . The energy lev- 
els are determined from the secular equation derived from 

HMD. 



Let us analyze now the role of the regularization op- 
erator, to remove the singularities in derivatives of the 
Green function as r,r' — > in ([26]). Performing the par- 
tial derivative with respect to x\ and taking first the limit 
r' — > yields 
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and £ = (E — Eo)/(hu>). The short-range singular be- 
havior of Hi (r) is given by the small-i asymptotics of the 
integrand in (|28|) . We calculate expansion of the inte- 
grand at t = 0, and later subtract the diverging terms 
giving rise to the singularity of i?;(r) at r — * 0. This 
allows us to decompose i?;(r) into a regular part Fi(t) 
and the rest which is singular at r — > 



Hi(t)=Fi(t)- 



where 
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With this result in mind one can now easily calculate Aki ■ 
Calculating partial derivative of (|27p with respect to Xk, 
and later applying the regularization operator, we ob- 
serve that the only nonvanishing terms comes from i*z(r), 
whereas the singular terms in Eq. (|28p are removed by 
d 3 r 3 . This demonstrates that the regularization opera- 
tor requires inclusion of the partial derivative of the form 
g r -2i+i in order to remove all the singular terms resulting 
from derivatives of the Green function. In this way we 



arrive at 

Aki (E) = -^-^H^k Fl (E,0)5 lk (31) 

where -Fi(O) is well defined, since the regularized inte- 
gral (I30j) now converges in both integration limits when 
r = . Substituting (|3"Tj) into set of coupled equations 
(|25p determining the eigenenergies, we observe that they 
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decouple, reducing to 



VxVyVz 



F k (E,0). 



(32) 



The energy spectrum generated by Eq. (|3"2"|) contains of 
three independent branches, that can labeled by k — 
Xjjj, z. This differs from the s-wave interacting atoms 
[40j |. where only one such a branch is found. We note 
that different branches do not describe to the excitations 
in the direction k, but rather results from the degeneracy 
of the p-wave states with respect to the projection of the 
angular momentum on the axis of quantization. We point 
out, that the total energy spectrum contains, apart from 
the eigenenergies determined from Eq. (|32p. also all the 
remaining eigenenergies of the noninteracting harmonic 
oscillator that have vanishing p-wave component at small 
r, and as a consequence are not affected by the p-wave 
pseudopotential. 

In the following we focus on the case of an axially sym- 
metric trap with lo x = oj v = lo+ = tjuj and uj z = u>. 
In such a geometry the projection of the total angular 
momentum on the z-axis is conserved, and the different 
branches in (|32p can be labeled by the projection quan- 



tum number m. First, we consider the case when 77 ia an 
integer or an inverse of integer, when the integral defining 
-Fi(O) can be given in a closed form. 

For pancake-shape traps, with r\ = 1/n, the set of equa- 
tions (1521) reduces to 
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where £' = E/(hcu), d z — y/h/(fj,ui z ), T(x) denotes the 
Euler Gamma function, and m is the projection of the 
angular momentum on the z-axis. We note that in the 
presence of an axially symmetric trap the degeneracy be- 
tween states with m = and m — ±1 is lifted. In order to 
evaluate the integral (f3"0"l) for 77 = 1/n we have applied the 
following identity l/{l- e -~ x / n ) = YX=l e~ xfe/ ™ / \l-e~ x ) 
(c.f. Ref. El). 

The closed formulas can be also derived for cigar-shape 
traps with integer r; = n: 



r,3 



r(j-f) 6 r (|-f) 
r(-i-f) r(i-f) 



"2 (n 



E 

fc=l 



} F 1, - 



^•1 

2 ' 



Si. -i 
2 i 



, 1 U 2 

> i TT~ 

r(-i-f) 



-6 



2 > 



n-1 



L U 2 J fc=i 



r 4 2 ' 4 2 ' 1 



m 



0. 



m = ±1, 



(35) 
(36) 



where F(a,b; c; x) denotes the hypergeometric function 
(4lj . As can be easily verified, for particular case of 
spherically symmetric trap (77 = 1), Eqs. (j53"|) - (f3"6")) co- 
incide with the result of Refs. [13,[2l[. To arrive at the 
Eqs. (f3"5| - (f3"6"|) we have applied a similar technique as in 
[4g|. making use of the following identities: 1/(1— e~ tn ) = 
kYl=l( l - e-*-*^*/")- 1 , and e' t( - n ^/(l - e"*") 2 = 

1 y-»n-l -i2nk/n _ „-t-i2itk / n\2 



directly from the integral (|3H)l 

F,(f;,0)-F z (£-2 J7> 0) = 
J F x (F,0)-2F ;c (£;-2r,,0) 



12/7 



r(l-f) 



r(-|) 

^(F- 477,0) 

_ 6r/ 2 r (f-f) 



r(i 



2 2 / 



(37) 



(38) 



The above equations together with (|3"0")) allows to evaluate 
F x .y,z(E,0) for arbitrary values of the energy and the 
trap anisotropy 77. 



In general, when r\ is not an integer nor inverse of an 
integer, Fi(E,0) can be calculated numerically from the 
integral representation (|30| for £ < 0, where the integral 
converges. This yields the part of spectrum representing 
the bound states of the interaction potential. In analogy 
to s-wave interactions 40], for energies £ > 0, above 
the zero point oscillation energy E one can apply the 
following recurrence relations that can be easily derived 



B. Energy spectrum for a model square- well 
interaction 

Results presented in the previous section have been de- 
rived using a zero-range potential, parameterized by the 
p-wave scattering volume a^. For s-wave interactions, in 
the regime of ultracold collisions, the scattering cross- 
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section acquires a constant value, and the pseudopoten- 
tial can be described by a single, energy-independent pa- 
rameter, e.g. s-wave scattering length. This however, 
does not hold when the interactions are dominated by 
the higher partial waves contributions, e.g. for polar- 
ized fermions, that can only interact via p-wave in the 
low-energy regime. In the case of p-wave interactions 
the range of applicability of the energy independent scat- 
tering length is very narrow and typically one has 
to apply the energy-dependent pseudopotential with the 
scattering length determined by the phase shifts of the 
short range potential. 

This behavior can be explained by analyzing the effec- 
tive range expansion of the corresponding p-wave phase 
shift: fc 3 cot(5i(fc) = -1/al - k 2 /(2R*) + C(fc 4 ), where 
R* can be interpreted as an effective range for p-wave 
interaction. The second term in the expansion can be 
neglected when fc|a p | (fci?*) 1 / 3 , which combined with 
the condition for the applicability of the pseudopotential: 
kR* <C 1, gives k\a p \ <C 1. Hence, the regime of applica- 
bility of the energy-independent p-wave pseudopotential 
is very limited, in contrast to the s-wave interactions, 
where the similar condition takes form fc|a s | <§C (fci?*) -1 , 
where a s is the s-wave scattering length and (kR*)^ 1 is 
large. 

In this section we consider a simple model of atoms 
interacting via square well potential. In the case of 
spherically symmetric trap the eigenenergies can be 
evaluated analytically, which allows to verify the accu- 
racy of the pseudopotential treatment. For the energy- 
dependent pseudopotential, we evaluate the energy levels 
self-consistently from Eqs. ([33]) - (f36|) . where a 3 is sub- 
stituted with the energy scattering volume a p (k) 3 = 
— tan 5 1 (k)/k 3 , where the p- wave phase shift 5i(k) is eval- 
uated for a square well potential. For negative energies 
this requires performing the analytical continuation of 
the phase shift, which can be obtained by substitution 
k = in , with E = -/i 2 k 2 /2/z [13]. 

Fig. [1] compares the exact energy levels for the square- 
well interaction, with the pseudopotential treatment 
assuming energy-dependent, and the standard energy- 
independent scattering volume. The comparison is per- 
formed for the two diameters of the square well: Ro — 
0.2d, and Rq = d, where d = y^h/fxoj. The depth U of the 
square well is varied in the vicinity of a resonance, that 
occurs when the bound state enters the potential well 
(c.f. the upper panels). We observe that EDP is very ac- 
curate, even when the size of the atom-atom interaction 
is comparable to the harmonic oscillator length deter- 
mining the characteristic scale of the external potential, 
and some deviations can be observed only for deep bound 
states in the case of Ro = d. On the other hand we note 
that the standard, energy-independent pseudopotential 
fails almost completely, and it works only very far from 
the resonance. 

Having established the regime of applicability of the 
pseudopotential method, we turn now to the case of 
anisotropic traps, where the exact levels are more difficult 
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FIG. 1: (Color online) Upper panels: The scattering volume 
dp versus depth U of the square well potential with diameter 
Ro = 0.2d (left panel) and Rq — d (right panel). Lower pan- 
els: energy spectrum for two fermions interacting via the same 
square-well potentials, and confined in isotropic harmonic 
trap of frequency u, versus square- well depth U. The exact 
results (black solid) are compared with the zero-range poten- 
tial model assuming energy-dependent (EDP, red dashed) and 
energy-independent (blue dotted) scattering length. Here, 
d — ^yh/JIuj denotes the harmonic oscillator length. 

to calculate, and we perform the analysis solely within 
the pseudopotential treatment. Figs. and [3] shows the 
energy levels calculated for 77 = 10 and 77 = 1/10. In con- 
trast to the previous analysis, here we plot the results as a 
function of the p-wave scattering volume, determining it 
for each depth U of the square- well at E = 0. In the case 
of energy-independent pseudopotentials this corresponds 
exactly to the energies given by Eqs ([55|) - ([5o| . respec- 
tively. This way of presenting the results let us to analyze 
the origin of the loop-like structures that can be observed 
in the part of the spectrum describing the bound states 
for calculations with the energy-independent pseudopo- 
tential. We observe that such structures are artificial, 
and disappear completely when one resorts to the cor- 
rect, energy-dependent pseudopotential. Similarly as in 
the case of isotropic trap, we observe that predictions 
of energy-dependent and energy-independent pseudopo- 
tentials differ almost completely, except the vicinity of 
points when a 3 = 0. 

C. Energy spectrum in the vicinity of a p-wave 
Feshbach resonance 

In this section we investigate the energy spectrum 
of two trapped spin-polarized fermions, with interac- 
tions modified by means of a p-wave magnetic Fesh- 
bach resonance. We start with deriving a formula de- 
scribing the energy-dependent p-wave scattering volume 
V(k) = — tanSi(k)/k 3 for a Feshbach resonance. To this 
end we utilize the multichannel quantum-defect theory 
(MQDT) [31], [13, HH, which has proven very effective 
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FIG. 2: (Color online) Energy spectrum for two fermions 
interacting in p-wave and confined in axially symmetric har- 
monic trap with r\ = uj±/uj z — 10, versus the scattering 
volume dp scaled by the harmonic oscillator length d z = 
y/h/JuJ^. The left and right panels show the eigenenergies for 
different values of the z-component of the angular momentum. 
The black-solid lines presents the spectrum calculated assum- 
ing the energy independent scattering volume, whereas the 
red-dotted lines are obtained for energy-dependent scattering 
volume cip(fc) of the square-well potential with Ro = 0.2d z 
(see text for details). 
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FIG. 3: (Color online) Same as Fig. [2] but for r] = 0.1. 



as few-parameter approach for describing scattering and 
bound states in neutral-neutr al l43l . liij , charged-neutral 
[45l |46|. and charged-charged |3l| particle collisions. 

In brief, MQDT is based on the separation of the short- 
range and the long-range behavior of the scattering wave 
function. The short-range wave function is insensitive 
to the total energy E, and can be described with a few 
quantum-defect parameters, that are insensitive both to 
the energy and the relative angular-momentum of the 
colliding particles [H, 53] • On the contrary, at large dis- 
tances the wave function vary rapidly with E, together 
with other scattering quantities defined at r — > 00 (e.g. 
phase shifts). 

In our treatment of a p-wave resonance, we consider for 
simplicity only one open and one closed channel. This is 
similar to CI (configuration interaction) model of a Fes- 
hbach resonance [48j, where a single closed channel rep- 
resents the effects of all the closed channels contributing 
to a resonance. In the two channel case, the application 
of MQDT is straightforward, and after some simple alge- 



bra (see for instance [H!]), one obtain the following phase 
shift S(E, I) in the open channel [50] 



6(E, I) — Sb s (E, I) — arctan 



\E,l) 



E - E n + L tan X(E, I) 



(39) 

The first term 6bg(E,l) is the phase shift resulting from 
the scattering in the open channel only, called the back- 
ground phase shift. The second term describes the res- 
onant contribution originating from a bound state in 
the closed channel, with energy E n located close to the 
threshold of the open channel. The width of the reso- 
nance T is multiplied by the MQDT function C~ 2 (E) 
that accounts for a proper threshold behavior as k — ► 0. 
In MQDT C~ 2 (E) relates the amplitudes of the solutions 
with short-range and long-range WKB-like normalization 
[H[53. The second MQDT function, tan X(E), describes 
a mixing between the two linearly independent solutions 
that occur in the quantum regime when we transform 
from the short-range to the long-range normalized wave 
functions [33|, [El]]. When k — > and the semiclassical 
approximation become valid at all distances, we have 
tanA(-E) -> and C~ 2 (E) -> 1, and the result ([39]) re- 
duces to the Fano-Beutler formula for the phase shift near 
a resonance [52j . 

We can assume that the energy of a bound state, varies 
approximately linearly with the magnetic field B: 



E n (B)=5fi(B-B n ), 



(40) 



where B n is the magnetic field when the bound state 
crosses the threshold of the open channel, and S/i is the 
difference of magnetic moments of the open and closed 
channels. In the case of a power law interaction potential 
r~ s , the phase shift and MQDT functions exhibit the 
following Wigner threshold behavior [53| 



S hg (E) t$ A,h 



21+1 



C~ 2 {E) E -^ B t k 2l+1 
tanA(£) ^°tanA(0) 



21 + 1 < s - 2 (41) 
all 1 (42) 
all 1 (43) 

The threshold behavior given by (|4"Tj) - (|4"3"l) is determined 
by the long-range part of the interaction potential, and 
for van der Waals forces the higher order terms in k can 
be neglected when energy is smaller than some character- 
istic energy i? v dw associated with C^/r 6 potential. For 
alkali -E v dw ranges from O.lmK to 30mK hence in 
the ultracold regime [E < 1/iK), one can safely use the 
approximations (|4l ]) -([l3" )) . Making use of Eqs. ([4H - ([43l) 
we can introduce the width of magnetic Feshbach reso- 
nance AB: 



km T - g " 2 (f> 
E->-0 2 ta,nri b g(E) 



-6fi AB, 



(44) 



and the resonance position Bo, that is shifted from B n 
due to the coupling between the open and closed channel 

r 

B = B n + — lim tanA(£). (45) 
25 fi -E->o 
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Now, substituting Eqs. 0DJ), gU), (gSJ) into ((39] 
performing some simple algebra we arrive at 

ta,n 5(E) = tan Sb g (E) 



and 



1 - 



AB [l + tan 2 Sb g (E)] 
B -Bo-fi + AB tan 2 6 hg (E) 



(46) 



For s-wave scattering, formula f4"6"j) yields the expres- 
sion for the energy-dependent scattering length (Eq.(75) 
in Ref. that apart from the scattering in very 

tight trapping potentials, reduces to the standard ex- 
pression a(B) = db g [1 — AB/(B — -Bo)]. In contrast, 
for p-wave resonances it is important to take into ac- 
count the energy-dependence of the resonance position 
and width. In this case the energy-dependent scattering 
volume V(E) reads 



V(E) = V hg (E) 



AB 



1 - 



B-Bo-^ + AB 



E 3 



(47) 



Here, E^ g = ti 2 /(2petb g ) is the energy associated with 
background scattering length db g = lim^^o Vbg(-E') 1 / 3 . 
Typically Et> g ~ E v dw, except the case when the 
background scattering volume is large: Vb g 3> R v dw, 
compared with the van der Waals length i? v dw = 
i(2^C 6 /ft 2 ) 1 / 4 . When E hg ~ E vdw , and for E < E vdw , 
when one does not need to take into account the effective- 
range expansion of Vb g (E) , Eq. (|47| can be simplified to 



V{E) S V hg 



1 - 



AB 



B - B - E/Sfx 



(48) 



We observe that in comparison to s-wave resonances, the 
position of the resonance is shifted for finite energies of 
colliding particles, and as we show later, this shift is sig- 
nificant for atoms confined in the tight traps (oj > 1kHz). 

Application of the result ([48]) in the Eqs. (|25|) and 
(I26[) yields the energy spectrum of the two spin-polarized 
fermions versus magnetic field. We have calculated the 
energy levels for a specific example of 40 K atoms, where 
p-wave Feshbach resonances has been observed experi- 
mentally and studied theoretically. In Ref. [HI Feshbach 
resonances have been parameterized in the sprit of an 
effective-range expansion 



k 3 cot 5x(k,B) = - 



1 



v(B) R(B) 



(49) 



The zero-energy scattering volume v, and the effective 
range R were approximated by the second-order polyno- 
mials in B 



1 



v(B) 
1 

RjB) 



c^B + c^B 2 



= d^ +d^B + d^B 2 



(50) 
(51) 



\mi\ 


c (0) (a.u.) 
dW(a.u.) 


c W(a.u.) 
d (1) (a.u.) 


c (2) (a.u.) 
d( 2 »(a'.u.) 





8.68155 x 10~ 5 
-1.64805 


-8.29778 x 10~ 7 
0.01523 


1.97732 x 10~ 9 
-3.54471 x 10~ 5 


1 


7.83424 x 10 -5 
-2.36792 


-7.45662 x 10~ 7 
0.02264 


1.76807 x 10" 9 
-5.45051 x 10~ 5 



TABLE I: Expansion coefficients of Eqs. ([50]) and ([ST]) 
for a p-wave Feshbach resonance between two 40 K atoms 
in / = | , mf — — |) state at Bo ~ 198G determined in 
Ref. [54| . The coefficients depend on the projection quan- 
tum number mi of the relative angular momentum along the 
axis of magnetic field. 



|mi| 


Bo[G] 


AB[G] 


Sfi/h [kHz/G] 


ab 9 [ao] 





198.85 


-20.342 


93.093 


-104.26 


1 


198.37 


-22.956 


92.667 


-99.539 



TABLE II: Resonance position Bq, resonance width AB, the 
difference of magnetic moments Sfi and the background scat- 
tering length aibg = (Vb g ) 1//3 , for a p-wave Feshbach resonance 
between two 40 K atoms in |/ = |,m/ = — |) state, depend- 
ing on their relative-angular-momcntum projection quantum 
number mi along the axis of magnetic field. The numbers 
have been extracted from the data in [5^] by fitting to the 
resonance formula (|48|l . 



where the expansion coefficient have been determined 
by fitting to the multichannel numerical calculations, re- 
maining with well agreement with the experimental mea- 
surements |54j . The expansion coefficients, listed in Ta- 
ble [T] depend on the projection of the relative angular 
momentum, mi, on the axis of the magnetic field. This 
multiplet structure of a p-wave Feshbach resonance re- 
sults from the small anisotropic spin dipole-dipole inter- 
action, that couples different partial waves, and breaks 
the rotational invariance. 

We have fitted expansion (|49]) to our formula , in the 
energy range < E < 500kHz and magnetic field range 
195G < B < 205G, obtaining the resonance parameters 
listed in Table HT] The agreement between (gSJ) and l[4"g]) 
remained on the level of few percent, which was similar 
to the accuracy of the effective range formula d4"9]) , when 
fitting to the numerical calculations in Ref. [54( . We have 
also tried the fit to the more accurate expression (I47p . 
including the effective range expansion for Vb s (E), but 
this has not further improved the quality of the fit. 

Figures [4] and [5] presents the energy spectrum for two 
40 K atoms, in the vicinity of a p-wave resonance at 
Bq w 198G. For simplicity we have chosen direction 
of applied magnetic field along the 2:-axis, in order to 
preserve the axial symmetry of the system. In such ge- 
ometry, the eigenstates with different mi are not coupled, 
and the resonance occurs at different value of magnetic 
field for mi = and \mi\ — 1, respectively. The reso- 
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B[G] 

FIG. 4: (Color online) Energy spectrum of two interacting 
fermions in the vicinity of a p-wave Feshbach resonance, ver- 
sus magnetic field B aligned along z direction. The atoms are 
confined in an axially-symmetric trap with r\ = u±_ /lu z — 1/5 
and uiz = 50 kHz. The blue-solid and red-dashed lines corre- 
spond to different projections of the angular momentum along 
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FIG. 5: Same as Fig. Q but for a cigar-shape trap with 
r\ — uij_/lu z = 5. 



nance can be observed as avoided crossing between the 
trap states represented by the horizontal lines, and the 
bound state corresponding to the skew line in the spec- 
trum. The effects of the trap anisotropy can be observed, 
as the different width of avoided crossing, depending on 
the excitation quantum number of the trap state (c.f. 
the avoided crossings with the 5th and 6th trap excited 
state for \mi\ = 1 in Figure [5]). We note that similarly 
to the square-well model, the inclusion of the energy- 
dependent scattering volume, eliminates the presence of 
the loop- like structures in the energy spectrum below the 
zero-point-oscillations energy. 



IV. SUMMARY AND CONCLUSIONS 

Summarizing, we have presented analytical solutions 
for two spin-polarized fermions interacting in p-wave and 



confined in a harmonic potential of arbitrary symmetry. 
For certain class of axially symmetric traps, with the as- 
pect ratio r\ integer or inverse of an integer, we have de- 
rived closed formulas for the energy spectrum. In gen- 
eral the eigenenergies can be easily evaluated numerically 
from some simple integral for energies smaller than the 
zero-point-oscillations energy: E < Eq, or with the help 
of recurrence relation for E > Eg. 

As an intermediate step we have discussed in details 
derivation of a generalized pseudopotential for all partial 
waves, introduced in Ref. [23j |. In the present paper we 
use slightly different form of the p-wave zero-range poten- 
tial in Cartesian coordinates, that turned out to be more 
convenient for our analytical calculations. In addition, 
we have presented the way to generate similar pseudopo- 
tentials for higher partial waves, showing the result also 
for the d-wave pseudopotential. Because of several dif- 
ferent formulations of the zero-range potentials for £ > 
existing in the literature, we have discussed how our pseu- 
dopotential connects to other approaches, and we have 
specified the conditions when different pseudopotcntials 
are equivalent. 

The higher partial wave interactions are vastly differ- 
ent from the s-wave one, and in the ultracold regime 
they typically require application of an energy-dependent 
pseudopotential (EDP). We illustrated this first for a 
model square-well potential, comparing the pseudopo- 
tential method with exact energy levels. We observed 
that EDP works well even when the range of atom-atom 
interaction is comparable to the characteristic size of an 
external potentials, which is well beyond the applicabil- 
ity of an ordinary pseudopotential. Then, we considered 
interactions of a pair of atoms close to a p-wave Fesh- 
bach resonance. In the framework of the multichannel 
quantum-defect theory, we have derived relatively sim- 
ple formula for the scattering volume in the vicinity of 
a resonance, that is quite similar to its s-wave analogue, 
which is well known in the literature. In this way we 
obtained the energy spectrum of two trapped fermions 
near a p-wave Feshbach resonance, that contrary to pre- 
dictions based on the energy-independent pseudopoten- 
tial, does not contain characteristic loop structures in 
the spectrum. In fact, the magnetic field dependence of 
the spectrum is quite similar to the one for s-wave in- 
teractions, in the presence of tight confining potentials 

We stress that our approach allows in principle for de- 
termination of eigenstates, however the resulting wave 
functions would exhibit the short-range 1/r 3 divergent 
behavior, and thus they could not be normalized in the 
standard way. The remedy is the application a delta-shell 
potential [2fj], that introduces a short-range cut-off that 
can be associated with the range of the atom-atom inter- 
action. Another solution is based on the regularization 
of the scalar product, according to the theory of gener- 
alized functions, or equivalently introducing a modified 
scalar product to define the normalized wave-functions 
[55l | . Moreover, for narrow p-wave Feshbach resonances 
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the contribution of the closed channel to the total wave 
function is significant, and cannot be neglected in the 
normalization. We stress that MQDT allows for rigor- 
ous determination of closed-channel component, which 
will not affect the results presented here for the spec- 
trum, nevertheless should be taken into account while 
determining eigenstates. These extensions, however, are 
beyond the scope of the present paper focusing solely on 
the properties of the energy spectrum, and will be inves- 
tigated elsewhere. 
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that diverge in the way described in Eq. (|A2 
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By comparison with (|A2jl . the finite part in the sense of 
Hadamard reads 



APPENDIX A: HADAMARD FINITE PART 
REGULARIZATION 



In this section we present a reasoning that leads to a 
finite part regularization of Eq. We basically fol- 

low the approach presented in the work of Estrada and 
Kanwal [3J]. Let's consider the following integral that 
diverges at the lower limit of integration 



Pf / dx 



ct>{x) 



dx 



fc-2 



4>(x) 



dx 



fc-l 



<M J) (Q) 3 



3=0 



(fc-l)! 

Substituting a — e and taking limit e — » + we obtain 



/oo 
dx 



4>(x) 



(Al) 



At small e, the integral exhibit the following behavior 



fc-i 



F(e) = F Q (e)+ coins + 

3=1 £ 



(A2) 



Pf 



1 



lim 

;=,0+ 



H{r-s) ^ (-lyS^jx) 

x k jUj - k + l)e fe -i- 1 

j=Q J VJ / 



(-l^-HneS^-^ix) 



(fc-l)! 

Here, H(x) denotes the Heaviside function 



H(x) 



x < 

1 x > 



(A6) 



(A7) 



with some constants c,-, and Fq(e) that is regular in the 
limit e —* 0. The finite part in the sense of Hadamard is 
defined as (see e.g. [35}) 



APPENDIX B: EXPRESSING THE 
PSEUDOPOTENTIAL IN TERMS OF Y lm (d) 
OPERATOR 



Pf QT dx = lim F (s). (A3) 

Splitting the integration into two intervals: (e,a) and 
(a, oo) with some arbitrary a > e, we can extract terms 



In order to prove identity (fT"5|) we calculate the action 
of the distribution T — 8^(r)Yi m (r)/r 2 on a test func- 
tion <f> that is regular and differentiable to arbitrary order 

/>oo 

(4>,T)= / drS^(r)f lm (r), (Bl) 
Jo 
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where 



fi m (r)= dnY lm {r)<f>(r). 



(B2) 



Performing Taylor expansion of cf>(r) around r = 0, one 
can easily obtain the small-r behavior of fi m (r) 



fi m (r)=ci m r l +0{r t+i ), r^O 



(B3) 



with 



(-1, 



i JdQY lm (r) ((fV)^) (0) (B4) 



In the next step we substitute (|B3|) into Eq. (jBip and 
make use of the property (fT4|) 

(0,T) - {-l)H\c lm 
4tt(-1)'/! 



Yi m (d)^2ci> m >r l 'Yi> m (r) 



{2l + l)\\ 



r=0 

(B5) 



This finally yields 



4ttZ! 



T) = — — (0, y Iro (a)5(r)) , (B6) 
which is equivalent to (| 1 5[) . 



APPENDIX C: SUMMATION FORMULA FOR 
Yim(d) OPERATORS 



The proof of Eq. (|19|) can be conveniently conducted 
using the Fourier transformation. We start with 



£ r ;m (d)Y?„(d')/(r)5(r') = 

i=-i 



,— i(kr+qr') 



m= — l 



(Cl) 



where /(r) and <7(r') are some test functions with Fourier 
transforms 



/» = 



■e lkr /(r), 



(2tt) 3 



(C2) 
(C3) 



The differential operator acting on the exponential func- 
tion produces: lz m (<9)e~' ikr = (— i) l Yi m (k)k l . Substitut- 
ing this result into (|C1[) and making use of the standard 
summation formula for the spherical harmonics (see for 
instance (4l| ) 



^m(k)^(q)=P/(kq), 



(C4) 



we arrive at 



Y lm (d)Y^(d')f(T) 9 (T') = 

in— — I 

2Z + 1 /■ fd 3 kd 3 q n/ * ,, 7j j 



4tt(-1)' 7 7 (2tt) 6 



fl(kq)fcV/(k)ff(q)c 



— -i(kr+qr ) 



(C5) 



Substituting the explicit form of the Legendre polynomial 
^U 35 ) — Si!=o c-kX l ~ 2k , with Cfc defined in ([20)) . we obtain 



m=—l 



Y lm {d)YL{d')f{r)9{v') = 

[A] 

= ^ c fc (V • V')'" 2 * V 2fc (V') 2fe /(r)ff(r'). (C6) 



fc=0 
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